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0\ '. ABSTRACT 

We have simulated the formation of an X-ray cluster in a cold dark matter universe using 12 different codes. 
The codes span the range of numerical techniques and implementations currently in use, including SPH and grid 
methods with fixed, deformable or multilevel meshes. The goal of this comparison is to assess the reliability of 
cosmological gas dynamical simulations of clusters in the simplest astrophysically relevant case, that in which the 
gas is assumed to be non-radiative. We compare images of the cluster at different epochs, global properties such as 
mass, temperature and X-ray luminosity, and radial profiles of various dynamical and thermodynamical quantities. 
On the whole, the agreement among the various simulations is gratifying although a number of discrepancies exist. 
Agreement is best for properties of the dark matter and worst for the total X-ray luminosity. Even in this case, 
simulations that adequately resolve the core radius of the gas distribution predict total X-ray luminosities that 
(^ • agree to within a factor of two. Other quantities are reproduced to much higher accuracy. For example, the 

\^ I temperature and gas mass fraction within the virial radius agree to about 10%, and the ratio of specific kinetic 

to thermal energies of the gas agree to about 5%. Various factors contribute to the spread in calculated cluster 
VO ' properties, including differences in the internal timing of the simulations. Based on the overall consistency of 

^D I results, we discuss a number of general properties of the cluster we have modelled. 

Subject headings: cosmology: theory — dark matter — galaxies: clusters — large-scale structure of universe 
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1. INTRODUCTION N-body techniques to follow the clustering evolution of a dissi- 
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cosmology. Two decades after the first cosmological Simula 
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that the validity of analytic approximations is often gauged by 
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reference to simulation results. 

The main limitation of N-body techniques is, of course, that 
they are relevant only to the evolution of dark matter In order 
to model the visible universe, it is necessary to include addi- 
tional physical processes. First, it is necessary to model a gas 
component that is gravitationally coupled to the dark matter In 
the simplest case, the gas may be assumed to be non-radiative. 
Although still highly simplified, this case has immediate appli- 
cations in the study of the hot plasma in galaxy clusters. At the 
next level of complexity, heating and cooling processes must 
be included. This is required, for example, to investigate the 
physical properties of the gas clouds responsible for QSO ab- 
sorption lines. Additional processes, such as star formation and 
the associated feedback of energy and mass, are necessary to 
model galaxy formation. 

Since the late 1980s a variety of techniques have been de- 
veloped to simulate gas dynamics and related processes in a 
cosmological context. In part inspired by the success of the N- 
body program, the first gas dynamical techniques were based on 
a particle representation of Lagrangian gas elements using the 
Smooth Particle Hydrodynamics (SPH) technique (Lucy 1977, 
Gingold & Monaghan 1977, Evrard 1988). Soon thereafter, 
fixed-mesh Eulerian methods were adapted (Cen et al. 1990, 
Cen 1992) and, more recently, Eulerian methods with submesh- 
ing (Bryan & Norman 1995) or deformable moving meshes 
(Gnedin 1995, Pen 1995, 1998) have been developed, as well as 
extensions of the SPH technique (Shapiro et al. 1996). These 
codes are actively being applied to a variety of cosmological 
problems, ranging from the formation of individual galaxies 
and galaxy clusters to the evolution of Lyman-a forest clouds 
and the large-scale galaxy distribution. 

Because of the inherent complexity of gas dynamics in a cos- 
mological context, such simulations are more difficult to vali- 
date than N-body simulations. Standard test cases with known 
analytical solutions (such as shock tubes) are far removed from 
the conditions prevailing in cosmological situations where the 
gas is coupled to dark matter and this, in turn, evolves through 
a hierarchy of mergers. The closest analogue to a realistic cos- 
mological problem is Bertschinger's (1985) solution for the col- 
lapse of a spherical cluster. Although this model provides a 
useful test of numerical hydrodynamics implementations, it ig- 
nores the merging processes that are a dominant aspect of the 
formation of reahstic clusters. In general, the strongly non- 
linear and asymmetric nature of gravitational evolution in a cos- 
mological context differs greatly from the regime that can be 
studied analytically or in laboratory experiments. 

In this paper we carry out an exercise intended as a step to- 
wards assessing the reliability of current numerical studies of 
cosmological gas dynamics. We address one of the simplest as- 
trophysically relevant problems, the formation of a large cluster 
in a hierarchical cold dark matter (CDM) model, using a variety 
of codes that span the entire range of numerical techniques in 
use today. The cluster problem is relatively simple because, ex- 
cept in the inner parts, the cooling time of the gas exceeds the 
age of the Universe and so, to a good approximation, the gas 
may be treated as non-radiative. 

The aim of this exercise is to assess the extent to which ex- 
isting modelling techniques give consistent and reproducible 
results in a realistic astrophysical application. Our compari- 
son is, by design, quite general. We simply specify the initial 
conditions for the formation of a cluster and let different simu- 
lators approach the problem in the manner they regard as most 
appropriate. Our comparison therefore encompasses not only 



the hydrodynamics simulation techniques themselves, but also 
individual choices of boundary conditions, resolution, internal 
variables such as the integration timesteps, and even the defini- 
tion of cluster center. We are therefore able to address issues 
such as the reproducibility of the X-ray luminosity and surface 
brightness of simulated clusters. It is not our intention to test the 
accuracy of any individual code: all the codes used in this pa- 
per have already been extensively tested against known analytic 
solutions. An earlier comparison of a subset of the techniques 
considered here was presented by Kang et al. (1994b). These 
authors simulated a large cosmological volume and focussed on 
statistical properties of the large-scale structure, rather than on 
the non-linear properties of an individual cluster that concern 
us here. 

This project was initiated as part of the activities of the 
program on "Cosmic radiation backgrounds and the formation 
of galaxies" which took place at the Institute for Theoretical 
Physics in Santa Barbara in 1995. Most active groups in the 
field of cosmological hydrodynamics simulations agreed to par- 
ticipate. Initial conditions for the cluster simulation were set up 
as described in Section 2. The codes used in the comparison 
are briefly described in Section 3. Participants were asked to 
analyze their results with a suite of predefined diagnostics, in- 
cluding images, global properties such as mass and X-ray lumi- 
nosity, and radial profiles of the dark matter density, gas density, 
temperature, etc. The images and a comparison of quantitative 
results are given in Section 4. Our paper concludes in Section 5 
with a summary and discussion of results, including some gen- 
eral conclusions regarding the properties of the simulated clus- 
ter. 

2. THE SIMULATION 

We simulated the formation of a galaxy cluster in a flat CDM 
universe. The initial fluctuation spectrum was taken to have an 
asymptotic spectral index, « = 1, and shape parameter, F = 0.25, 
the value suggested by observations of large-scale structure (eg. 
Efstathiou, Bond & White 1992). The cosmological param- 
eters assumed were: mean density, 17 = 1; Hubble constant, 
Hq = 50km s"' Mpc"' ; present-day linear rms mass fluctuations 
in spherical top hat spheres of radius 16 Mpc, erg = 0.9; and 
baryon density (in units of the critical density), il;, = 0. 1 . 

2.1. Initial conditions 

Initial conditions were laid down using the constrained Gaus- 
sian random field algorithm of Hoffman & Ribak (1991). The 
cluster perturbation was chosen to correspond to a 3(7 peak 
of the density field smoothed with a Gaussian filter of radius 
ro = 10 Mpc [in exp(-0.5(r/ro)^)]. The perturbation was cen- 
tered on a cubic region of side L = 64 Mpc. We used the fit to 
the CDM transfer function given by eqn (G3) of Bardeen et al. 
(1986) and recommended a starting epoch of z = 20. 

To offer flexibility, the initial conditions were generated at 
very high resolution. Two alternative forms were supplied: 

(i) The dimensionless linear 5p/ p field (normalized to the 
present), tabulated on a 256"^ cubic mesh. 

(ii) The linear theory displacements for 256^ points on a cu- 
bic mesh. 

The initial conditions were generated by Shaun Cole. 
These are publically availab le on the Internet at h ttp://star- 
www.dur.ac. uk~csf/clusdata/ or by request from CSF. 
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2.2. Simulation diagnostics 

The simulation data were output at redshifts z =8, 4, 2, 1, 0.5, 
and and the following properties were calculated: 

Images. 

At each epoch, 2-D images were generated of various quanti- 
ties in the central (32 Mpc)^ comoving volume. The quantity of 
interest was projected along the z-axis, smoothed as specified 
below, and tabulated on a 1024^ mesh. Two smoothings were 
employed: 

(i) A Gaussian smoothing with kernel, exp(-0.5(r/ro)^), at a 
fixed resolution of ro = 250 kpc (comoving). 

(ii) A smoothing of choice, as determined by each simulator. 
The first set of images was used for a uniform comparison of 
all the models, while the second set was supplied in order to 
display the results of each technique in the best possible light. 

Images of the following quantities were made: 

(i) Projected dark matter density (in MoMpc"^) 

(ii) Projected gas density (in MqMpc"^) 

(iii) X-ray surface brightness, f Cxdl (X-ray emissivities 
per unit volume were calculated as Cx = p^T^I^, with p in 
MoMpc"^ and T in K.) 

(iii) Emission-weighted temperature {j CxTdl/ j Cxdl). 

Global properties. 

We defined the cluster to be the mass contained inside a 
sphere of radius, r2oo, such that the mean interior overdensity 
is 200. The following global properties of the cluster were then 
computed at z = 0: 
(i) The value of r2oo 

(ii) Mdm : total dark matter mass (in Mq) 
(iii) Vdm : rms velocity of dark matter particles (in km s"') 
(iv) Mgas : total gas mass (in Mq) 
(v) T : mean (mass weighted) temperature (in K) 
(vi) U : total bulk kinetic energy of the gas (in ergs) 
(vii) Ltot = Jq'"" CxdV {Cx units as above; V in Mpc-') 
(viii) I = ^j niiXiXi/ ^^ m,: inertia tensor for dark matter and 
gas (in Mpc^). 

Radial profiles. 

In 15 spherical shells of logarithmic width 0.2 dex and in- 
ner radii lOkpc < r < lOMpc, the following quantities were 
obtained: 

(i) Pdmif) '■ dark matter density profile (in MqMpc"^) 
(ii) o'dm(r) '. dark matter velocity dispersion profile (in km s"^) 
(ii) Pgas('") : gas density profile (in MqMpc"-') 
(iv) T(r) : mass-weighted gas temperature profile (in K) 
(v) Cxir) : "X-ray luminosity" profile calculated as the to- 
tal luminosity in each bin, divided by its volume (in units as 
above). 

Each simulator was provided with the two initial conditions 
files and the list of required diagnostics (in a prespecified for- 
mat). Everything else was left to the discretion of each simula- 
tor including, for example, the definition of the cluster center. 
All data were sent directly to the organizers (CSF and SDMW), 
and participants were strongly discouraged from private inter- 
comparison of results. A surprising number of iterations was 
required to obtain consistent outputs in a single set of units and 
formats. 

Wadsley joined the project after the original deadline had ex- 
pired and the first set of results was known. Discrepancies in a 
preliminary comparison of results led Gnedin to revise his code 



and resubmit a new simulation. Bryan's stated spatial resolu- 
tion was changed from 80 to 30 kpc after preliminary compar- 
ison suggested that he had been too pessimistic in stating his 
resolution. 

3. THE CODES 

The numerical codes used for this project employ a variety 
of techniques to solve the evolution equations for a two com- 
ponent fluid of dark matter and non-radiative gas coupled by 
gravity. In what follows, each code is identified by the name 
of the author who was primarily responsible for carrying out 
each simulation. The codes are of two general types: SPH and 
grid-based. SPH simulations were carried out by Couchman, 
Evrard, Jenkins, Navarro, Owen, Steinmetz and Wadsley. The 
simulations by Couchman and Jenkins used the same basic code 
(HYDRA), the serial version in the former case and a parallel 
version in the latter (These two simulations were done inde- 
pendently and used different numbers of particles and different 
values for the simulation parameters: gravitational softening, 
smoothing length, timestep, etc.) Owen's code differs from the 
others in the use of an anisotropic SPH kernel. The grid-based 
methods employ either a single, fixed mesh (Cen, Yepes), a 2- 
level multi-mesh (Bryan) or a deformable mesh (Gnedin, Pen). 
Warren carried out a high resolution simulation of the evolution 
of the dark matter only. 

A brief description of each code follows, together with ref- 
erences where the reader may find a fuller discussion of tech- 
niques and the tests to which each code has been subjected. 
Details of each simulation are given in Table 1 . 

3.1. SPH codes 
3.1.1. Couchman & Thomas - Hydra (Adaptive P^M-SPH) 

Hydra is functionally equivalent to the standard particle- 
particle-particle-mesh, N-body-SPH (P^M-SPH) implementa- 
tion, but with the automatic placement of a hierarchy of refined 
meshes in regions of high particle density. This avoids the dra- 
matic performance degradation caused by the direct summation 
(PP) component of standard P^M codes under heavy particle 
clustering. In the present simulation, at a redshift of z = 0.5, 
the cpu time per step had increased by a factor of 4.5 from the 
essentially uniform initial conditions at z = 49, and remained 
at this level to the end of the simulation. A maximum of four 
levels of mesh refinment was chosen by the code. 

The code automatically chooses a global timestep to ensure 
accurate time integration. This value is determined by the max- 
imum instantaneous values of particle velocities and accelera- 
tions of both gas and dark matter particles. An optimal low- 
order integration scheme is used for advancing particle posi- 
tions and velocities. Full details of the code are available in 
Couchman, Thomas and Pearce (1995) and the source code 
may be found in Couchman, Pearce and Thomas (1996). 

The supplied initial displacement field was degraded to the 
resolution used, 64-^ dark matter particles, simply by sampling 
every fourth position in each dimension. This rather crude 
method of resampling, although simple, has the disadvantage of 
introducing noise into the perturbed particle distribution above 
the effective Nyquist frequency of the 64^ particles. Dark mat- 
ter particles were displaced from a uniformly spaced grid and 
gas particles were placed on top of the dark matter particles. 
Particle displacements were scaled to correspond to a start red- 
shift of 49. The center of the final cluster was identified with 
the density peak within r2oo- 
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3.1.2. Jenkins & Pearce - Parallel Hydra 

This simulation used a parallel version of Hydra, the code 
just described in 3.1.1, whose distinguishing feature is its abil- 
ity to place high resolution meshes recursively around clus- 
tered regions. The SPH calculation used an M4 spline ker- 
nel containing an average of 32 particles. Details of the se- 
rial version of this code may be found in Couchman, Thomas 
and Pearce (1995), while details of the parallel implementation 
may be found in Pearce et al. (1995) and Pearce & Couchman 
(1997). 

The simulation was carried out on a Cray-T3D. Initial con- 
ditions were laid down by perturbing 128^ particles distributed 
in a "glass" configuration. This was generated in the manner 
described by White (1996), ie. by evolving a Poissonian distri- 
bution of points, with the sign of gravity reversed, over many 
thousands of expansion factors. To optimize the resolution in 
the region of interest, the computational volume was divided 
into two parts, a high resolution spherical region containing 1 /4 
of the volume and centered on the location of the constrained 
peak, and a coarsely sampled exterior region. Dark matter and 
gas particles (initially coincident) were placed in the high res- 
olution region, and dark matter particles only in the exterior 
region. The coarse sampling was achieved by smoothing the 
distribution with a nearest grid point (NGP) assignment on a 
64-' mesh, so that, on average, each particle was 8 times more 
massive than particles in the high resolution region. This pro- 
cedure reduced the particle number from 2097152 dark matter 
particles to 1247217 of both species, about one million of which 
lay in the high resolution region. The initial particle positions 
were set up at z = 20 using a trilinear interpolation of the dis- 
placement field in the 8 mesh points surrounding each particle. 
Velocities were assigned from the Zel'dovich approximation. 
The center of the final cluster was defined to be the position of 
the particle with the lowest gravitational potential. 

Since this and Couchman's simulation were carried out with 
the same code, one in parallel and the other in serial mode, any 
differences in the results must be due to differences in the ini- 
tial conditions or the choice of integration parameters. We have 
checked that running Couchman's initial conditions in parallel 
mode does not alter his results in any significant way. 

3.1.3. Evrard-P^M-SPH 

The P^M-SPH code combines the P^M code of Efstathiou 
& Eastwood (1981) with an adaptive kernel SPH scheme, as 
described in Evrard (1988). The simulation for this study em- 
ployed a two-level mass hierarchy, with a high resolution (64^ 
effective) inner zone of both dark matter and gas surrounded 
by dark matter at low resolution (32^). The A^^ mesh data were 
generated by NGP subsampling of the original 256^ displace- 
ment field. The mapping of the high resolution zone was deter- 
mined by a low resolution (32^) N-body simulation; particles 
within a final density contrast of 6 centered on the group in this 
run define a Lagrangian mask used to generate the two-level 
initial conditions of the full run. Masked locations in the 32^ 
subsampled field were locally "exploded" to a factor 2 higher 
linear resolution, generating an effective 64^ resolution within 
the non-linear parts of the cluster. This procedure assures no 
contamination of low resolution particles within the cluster in 
the production run. The run used 30456 particles for each high 
resolution component (gas and dark matter) and 28961 particles 
at low resolution, with a 128^ Fourier mesh for the long-range 
gravity. The center of the final cluster was defined by the most 



bound dark matter particle. 

The number of interacting neighbors within the smoothing 
kernel controls the hydrodynamics resolution of the calcula- 
tion. This parameter was set so that approximately 168 par- 
ticles lie within a sphere of radius 2/z around any particle. As 
discussed by Owen and Villumsen below, the value of this pa- 
rameter varies considerably among experiments in the litera- 
ture. The value employed here is larger than "typical" values 
and reflects a desire to minimize the Poisson noise inherent in 
the kernel summations required for calculation of the density 
and pressure gradient terms. 

The computation was performed on a local HP workstation. 
The modest memory and CPU requirements of this calculation 
reflect its nature as closer to "everyday" than "state-of-the- 
art". It is representative of the type of runs used in ensembles 
to investigate statistical aspects of the cluster population (e.g. 
Mohr& Evrard 1997). 

3.1.4. Navarro - Grape+SPH 

The code used was the N-body/SPH code described by 
Navarro & White (1993), adapted to compute gravitational ac- 
celerations using a GRAPE-3 board. The neighbor lists needed 
for the SPH computations are also retrieved from the GRAPE 
and processed in the front-end workstation. The implementa- 
tion of these modifications is straightforward and very similar 
to that described by Steinmetz (1996), where the reader may 
find far more details. 

The initial conditions were realized by perturbing a cubic 
grid of particles with the displacement field made available with 
the initial conditions package. The system was divided in two 
zones, an inner cube of size 38 Mpc which was filled with 40^ 
dark matter and 40^ gas particles, surrounded by a sphere of 
diameter 64 Mpc. The region outside the inner cube was filled 
with ^ 5,000 low-resolution dark matter particles of radially 
increasing mass. Initially, gas and dark matter particles were 
placed on top of each other and were given the same veloci- 
ties, computed using the Zel'dovich approximation. The final 
cluster center was calculated using a concentric sphere method 
that isolates the highest density peak iteratively by computing 
the center of mass of a sphere and successively removing the 
outermost particle, until only about 100 particles are left. 

The simulation was run on the SPARClO/GRAPE-3 system 
at Edinburgh University. 

3.1.5. Steinmetz - GrapeSPH 

This simulation was performed using GrapeSPH (Steinmetz 
1996), a direct summation hybrid N-body/SPH code especially 
designed to take advantage of the hardware N-body integrator 
GRAPE (Sugimoto et al. 1990). It is highly adaptive in space 
and time through the use of individual particle timesteps and 
individual smoothing lengths. Details of the code such as the 
adaptive smoothing length or the multiple time stepping proce- 
dure, are presented in Steinmetz & Miiller (1993) and in Stein- 
metz (1996). 

The simulation used a multi-mass technique similar to that 
described by Porter (1985). Firstly, a low resolution (32^ par- 
ticles) P^M simulation of the full periodic volume was per- 
formed. The initial conditions were drawn from the distribu- 
tion supplied by averaging positions and velocities in cubes of 
8^ particles. At z = the cluster which formed near the center 
was identified and its virial radius, r2oo, determined. Particles 
within r2oo were marked and traced back to redshift z = 20. A 
sphere was then drawn containing all these particles - the high 
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resolution region. Particles within that sphere were replaced by 
the corresponding particles drawn from higher resolution initial 
conditions. Thus, the particle number in the high resolution re- 
gion was increased by factors of 8-64. Particles outside the high 
resolution region were combined into larger mass nodes us- 
ing the tree-pruning technique described in Porter (1985). The 
mass of a particle thus increases logarithmically with distance 
from the central sphere. Starting from these initial conditions, 
the full simulation was performed, with gas dynamics followed 
only in the high resolution region. The center of the final clus- 
ter was defined as the center of mass of the smallest (radius 
125 kpc) of a series of 7 concentric spheres of progressively 
decreasing radius. (This agreed to about 5% with the minimum 
of the gravitational potential.) 

In the hardware integrator GRAPE, the interparticle force is 
hardwired to obey a Plummer force-law. Thus, periodic bound- 
ary conditions cannot easily be realized (for a more recent de- 
velopment, see Huss, Jain & Steinmetz 1998). Because of the 
tree-pruning, however, the CPU time scales only logarithmi- 
cally with box size for a given numerical resolution. A typical 
application thus starts from a very large simulation sphere as- 
suming vacuum boundaries. Since the computational box sup- 
plied for this cluster simulation was relatively small, effects due 
to the finite box size cannot be excluded and this may also affect 
the comparison with grid based methods. In order to minimize 
the effects of finite box size and vacuum boundary conditions, 
the simulation strategy was slightly modified. Tree pruning was 
not applied to the original box, but to an enlarged box including 
the 26 neighboring periodic replicas of the original. Hence, vac- 
uum boundaries apply to a surface of radius r = 1.5 /box, rather 
than r = 0.5 /box- 

A variety of simulations with differing numerical resolution, 
particle numbers and size of the high resolution region were 
performed. Results from one simulation only have been in- 
cluded in this paper. This probably reflects the best compromise 
between resolution and computational cost. In this simulation, 
which took about 28 hours of CPU and 22 MBytes of mem- 
ory, ^ 15000 gas and dark matter particles ended up within the 
virial radius of the cluster at redshift z = 0, a resolution simi- 
lar to that achieved by Evrard, Navarro, Couchman and Wad- 
sely. The largest simulation carried out had the same number 
of gas particles but 8 times as many dark matter particles. This 
run consumed a total CPU time of 254 hours and required 45 
MBytes of memory. 

3.1.6. Wadsley&Bond-P^MG-SPH 

The Wadsley and Bond (1997) P^MG-SPH code used in this 
cluster comparison combines SPH for the hydrodynamics with 
an iterative multigrid scheme to solve for the non-periodic grav- 
itational potential with a particle-particle correction for subgrid 
forces. A recursive linked list is used to locate neighbor parti- 
cles for SPH. At each timestep, a multipole expansion is used to 
obtain the gravitational potential boundary conditions on a 128^ 
grid. The multigrid technique is quite competitive with Fast 
Fourier Transform methods in speed and can more efficiently 
treat non-periodic configurations, for which the P^MG-SPH 
code is designed. It is typically used to compute highly active 
inner regions at high resolution, with large scale tides treated 
using a sequence of progressively lower resolution spherical 
shells in the initial conditions. The force is augmented by a 
measured external tidal field evolved using linear theory. 

This cluster had a Gaussian filter scale for the peak which 
was too large for the periodic box size to treat the tidal environ- 



ment adequately {i.e., the cluster was apodized). To a spherical 
high resolution region of radius 25 Mpc, we could only add a 
single lower resolution shell extending to 32 Mpc in radius. The 
25 Mpc choice was based on the region the peak-patch theory 
(Bond and Myers 1996) suggests would have collapsed. The 
low resolution particles had 8 times the mass of the high reso- 
lution ones. The self consistent linearly-evolved external shear 
was also applied to the entire region. There were 74127 gas 
and 74127 dark particles used in the simulation. The initial 
256^ displacement field was sampled at every fourth lattice site 
to transfer onto the computational grid for the high resolution 
region. A similar transfer was done for the low resolution re- 
gion, with slight smoothing added. Couchman used the same 
one-in-four transfer method, probably accounting for the sim- 
ilarities with the Wadsley and Bond result, especially with re- 
gard to timing. Discrepancies may be due, at least in part, to 
his use of periodic boundary conditions. The center of the final 
cluster was taken to be the center of mass of the largest group 
in the simulation identified with a standard friends-of-friends 
group finder. 

The computation was run on a Dec- Alpha EV5 and required 
119 CPU hours and 100 Mb of memory. The current version 
of this code has a significantly accelerated particle-particle sec- 
tion, using tree techniques (the PP section of the old gravity 
solver slowed by a factor of 10 by z = 0). The same computa- 
tion now takes 33.6 CPU hours and 50 Mb of memory. 

3.1.7. Owen & Villumsen - Adaptive SPH 

This simulation was performed with a variant of the SPH 
method called Adaptive Smoothed Particle Hydrodynamics, or 
ASPH. ASPH generalizes the isotropic sampling of SPH by 
associating an individual, ellipsoidal interpolation kernel with 
each ASPH node, the size, shape, and orientation of which is 
evolved using the local deformation tensor dvi/dxj. The goal of 
the algorithm is to maintain a constant number of particles per 
smoothing length in all directions at all times for each ASPH 
particle. This anisotropic sampling allows the ASPH resolution 
scale to better adapt to the local flow of material as compared 
with the isotropic sampling of traditional SPH, thereby max- 
imizing the resulting spatial resolution for a given number of 
particles. Another way of stating this is to say that in the frame 
defined by the kernel, the distribution is locally isotropic. 

The ASPH formalism is meaningful only when particles are 
treated not as particles but as moving centers of interpolation. 
The effects of momentum non-conservation are minimized so 
long as internal consistency in the ASPH kernel field is main- 
tained, a condition we enforce by renormalizing the ASPH ker- 
nels periodically. The prescription for this renormalization, the 
ASPH algorithm, and the code used here are described in detail 
in Owen et al. (1997), and an earlier discussion of ASPH may 
be found in Shapiro et al. (1996). 

The ASPH interpolation between nodes is performed using 
the bi-cubic spline interpolation kernel, which formally ex- 
tends for two smoothing scales, h. Beyond this point, the bi- 
cubic spline falls to zero, and therefore only nodes within 2h 
can interact with each other The local smoothing scales were 
initialized such that there are roughly 2 nodes per smoothing 
scale, so each node "sees" a radius of 4 nodes, or a total of 
4/3 TT 4-^ w 268 nodes. Of course, since the bi-cubic spline 
weighting falls to zero near the edge of this sampling volume, 
each node effectively interacts with only about 4/3 tt 3^ « 1 13 
neighbors. While this represents a much larger number of 
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neighboring particles than most contemporary SPH implemen- 
tations (a more typical value is 1 particle per h, yielding 32 
particles in a volume of radius 2h) it appears that keeping 1 .5- 
2 particles per h yields more reliable results for a wide vari- 
ety of hydrodynamical test problems (see also Balsara 1995). 
Evrard's SPH simulation also used a a large number of inter- 
acting particles. The disadvantage of this choice, of course, is 
that it decreases the effective resolution. An additional feature 
of the code is that it uses a compact, higher order interpolation 
kernel for the artificial viscous interactions in an effort to more 
closely confine the effects of the artificial viscosity to shocked 
regions. The gravitational interactions are evaluated using a 
straight, single-level Particle-Mesh (PM) technique. 

This experiment was performed using 32^ dark matter and 
32-^ ASPH particles, and a 256-^ PM grid for the gravity. Since 
there are equal numbers of ASPH and dark matter particles, 
each dark matter particle is 9 times as massive as an ASPH par- 
ticle. The initial conditions were generated for z = 20 by per- 
turbing the 2 X 32^ ASPH and dark matter particles (initially 
exactly overlaying each other) from a cubical lattice with tri- 
linear interpolation based upon the supplied displacement field, 
and assigning velocities using the Zel'dovich approximation. 
The center of the final cluster was defined through an iterative 
approach, similar to that used by Steinmetz. In this case, the 
radius of each successive sphere was shrunk by a factor of 0.9 
and the iterations were stopped when the center of mass shifted 
by less than a small tolerance. 

It is worth commenting on the fact that with only 32^ ASPH 
and dark matter particles, this is by far the lowest resolution of 
the SPH experiments presented in this paper. This is primarily 
due to the fact that at the time this experiment was performed 
our 3-D ASPH code had only just been completed, and this was 
one of the first problems tackled with that (then) highly experi- 
mental code. Due to the lateness of our entry into this project, 
we only had time to confirm that the experiment appeared to 
have run successfully and submit that initial run. The current 
version of our 3-D ASPH code is competitive with other con- 
temporary cosmological hydrodynamical codes. Despite these 
limitations, though, it is interesting to compare the results of 
this simulation with the others in order to quantify how well 
the regions which are resolved match the results of the higher 
resolution models. 

3.2. Grid-based methods 

3.2. 1 . Bryan & Norman - SAMR 

A newly-developed, structured, adaptive mesh refinement 
(SAMR) code was used to perform this simulation. This 
method was designed to provide adaptive resolution while pre- 
serving the shock-capturing characteristics of an Eulerian hy- 
drodynamics scheme. The code identifies regions requiring 
higher resolution and places one or more finer sub-meshes over 
these areas in order to better resolve their dynamics. There is 
two-way communication between a grid and its 'child' meshes: 
boundary conditions go from coarse to fine, while the improved 
solution on the finer mesh is used to update the coarse 'parent' 
grid. The grid placement and movement is done automatically 
and dynamically, so interesting features can be followed at high 
resolution without interruption. Since sub-grids can have sub- 
sub-grids, this process is not limited to just two levels. The 
control algorithm for advancing the grid hierarchy is similar 
to that suggested by Berger & Colella (1989), and the equa- 
tions of hydrodynamics are solved on each grid with a version 



of the piecewise parabolic method (PPM) modified for cosmol- 
ogy (Bryan et al. 1995). Dark matter was modeled with par- 
ticles, and gravitational forces were computed via an adaptive 
particle-mesh scheme. Poisson's equation was solved on each 
mesh using the Fast Fourier Transform. 

The simulation was initialized at z = 30 with two grids al- 
ready in place. The first is the root grid covering the entire 64 
Mpc^ domain with 64^ cells. The second grid is also 64-^ cells 
but is only 32 Mpc on a side and is centered on the cluster, 
yielding an initial cell size of 500 kpc (the initial conditions 
were provided at higher resolution but were smoothed with a 
sharp k-space filter where appropriate). The refinement crite- 
ria was based on local density, so any cell with a baryon mass 
of 3.5 X 1O^M0 or more, was refined, but only within the cen- 
tral, high-resolution region. At z = 0.5, the hierarchy consisted 
of more than 300 grids spread out over 7 levels of refinement. 
Each level had twice the resolution of the one above, produc- 
ing, in very small regions, a cell size of 8 comoving kiloparsecs. 
The final cluster center adopted was the cell-center of the cell 
with the highest baryonic density. The simulation was carried 
out on four processors of an SGI Power Challenge at the Na- 
tional Center for Supercomputing Applications (NCSA). 

3.2.2. Cen,Bode,Xu&Ostriker-TVD 

This simulation employed a new shock-capturing Eulerian 
cosmological hydrodynamics code based on Harten's Total 
Variation Diminishing (TVD) scheme (Harten 1983), and de- 
scribed in Ryu et al. (1993). The original TVD scheme was im- 
proved by adding one additional variable (entropy) and its evo- 
lution equation to the conventional hydrodynamics equations. 
This improvement eliminates otherwise large artificial entropy 
generation in regions where the gas is not shocked. Details of 
this treatment can be found in Ryu et al. (1993). The code is 
able to capture a strong shock within 1 -2 cells and a sharp den- 
sity discontinuity within 3-5 cells. Poisson's equation is solved 
using the Fast Fourier Transform. The code is accurate in terms 
of global energy conservation to about 1 %, as gauged by the 
Layzer-Irvine equation. 

Initial conditions were laid down on a 512^ uniform grid, us- 
ing the particle positions, velocities and gas densities provided. 
The initial gas temperature was set to a low value, 1 13 K. The 
total number of particles was 256^ and the total number of fluid 
cells was 512^. The simulation started at z = 40. The final clus- 
ter center was taken to be the cell with the highest X-ray lumi- 
nosity. 

The simulation was run on an IBM SP2 at the Cornell Theory 
Center Sixty-four SP2 processors were used for the simulation 
for about 83 wallclock hours with 600 timesteps. The code 
is well parallelized on the SP2 (as well as on the Cray-T3E) 
and achieves an efficiency of about 50% on 64 SP2 processors 
(Bodeeffl/. 1996). 

3.2.3. Pen - Moving Mesh Hydrodynamics 

The simulations were all performed with an early version of 
the Moving Mesh Hydrodynamics and N-body code (MMH for 
short; Pen 1995, 1998). Its main features are a full curvilin- 
ear TVD hydro code with a curvilinear PM N-body code on a 
moving coordinate system. The full Euler equations are solved 
using characteristics in explicit flux-conservative form. The 
curvilinear coordinates used in the code are derivable from a 
gradient of the Cartesian coordinate system. If x' are the Carte- 
sian coordinates, the curvilinear coordinates are ^' =x' + d^p(j){£). 
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The transformation is completely specified by the single poten- 
tial field 0(^, 0- During the evolution any one constraint can be 
satisfied by the grid. In our case, we follow the mass field such 
that the mass per unit grid cell remains approximately constant. 
This gives all the dynamic range advantanges of SPH combined 
with the speed and high resolution of grid algorithms. For rea- 
sons of cpu economy (the computational cost increases linearly 
with compression factor), this run was constrained to compress 
by at most a factor of 10 in length, or a factor of 1000 in density. 

The potential deformation maintains a very regular grid 
structure in high density regions. The gravity and grid defor- 
mation equations are solved using a hierarchical multigrid al- 
gorithm for linear elliptic equations. These are solved in linear 
time, and are asymptotically faster than the FFT gravity solver. 
At the same time, adaptive dynamic resolution is achieved. The 
gravitational softening of 45 kpc listed in Table 1 applies to the 
central region of the cluster; the average softening is 450 kpc. 
The final cluster center was identified with the minimum in the 
gravitational potential. 

The algorithmic cost per particle per timestep is very small, 
^ 300 FLOP (floating point operations). The cost for the grid 
deformation, gravity and hydrodynamics adds to about 20k 
FLOP per grid cell per timestep. If memory is available, we 
always use 8 particles per grid cell, at which point the particles 
only account for a small portion of the computation time. This 
ensures that we are unlikely to encounter artifacts due to 2-body 
relaxation. 

The code runs in parallel on shared memory machines with- 
out load balancing problems. The simulation was carried out 
using a 128"' grid and 256"' particles. Currently each timestep 
takes 60 seconds on a 16 processor Origin 2000 at 195 Mhz. 
The whole run takes 1600 time steps or about one day. At the 
time the actual simulations were performed, the best available 
machine was a 75 Mhz R8k Power Challenge, where the run on 
8 processors took 60 hours. 

Initial conditions were specified on an initially uniform grid. 
The fluid perturbation variables were set up on the grid, and 
particles displaced using the Zel'dovich approximation. (At the 
time of writing, the code no longer uses the Zel'dovich approxi- 
mation, but instead varies the mass of each particle.) The initial 
grid can now be adjusted to resolve hierarchically any region of 
interest with arbitrary accuracy. We expect the new version to 
perform significantly better. 

3.2.4. Gnedin - Smooth Lagrangian Hydrodynamics 

In this code, the hydrodynamical evolution of the gas is fol- 
lowed using the Smooth Lagrangian Hydrodynamics or SLH 
method (Gnedin 1995), in which all physical quantities are de- 
fined in quasi-Lagrangian space, q'^, and Eulerian positions, x', 
are considered as dynamical variables. The imaginary mesh 
connecting Eulerian positions, x', thus moves with the fluid un- 
til one of eigenvalues of the deformation tensor, A^ = dx^ jdc^, 
becomes smaller than the predefined softening parameter. A*. 
Then in the direction corresponding to this eigenvalue the mesh 
gradually decelerates and progressively approaches (but never 
fully reaches) the locally stationary mesh, until (and if) the cor- 
responding eigenvalue of the deformation tensor begins to in- 
crease. This process of softening of the Lagrangian flow pre- 
vents severe mesh distortions which can cause the stability and 
accuracy of a purely Lagrangian code to deteriorate. The grav- 
itational force in the code is computed using the P^M method 
and is subject to the Gravitational Consistency Condition as de- 
scribed in Gnedin & Bertschinger (1996). 



Initial conditions were set up by sampling the supplied fields 
on a 64^ mesh and using the Zel'dovich approximation to ad- 
vance the dynamic variables to z = 20. The cluster center was 
defined using the DENMAX algorithm. 

3.2.5. Yepes, Khokhlov & Klypin - PM-FCT 

The code used for this simulation is a combination of an 
Eulerian hydrodynamical code based on the Flux-Corrected- 
Transport (FCT) technique (Boris 1971, Boris & Book 1973, 
1976) and a standard Particle-Mesh N-body code (Kates, Ko- 
tok and Klypin 1990). It uses the "low phase error algorithm" 
whereby phase errors in convection are reduced on the uniform 
grid to fourth order (Boris & Book 1976, Oran & Boris 1986). 
This algorithm is applied to the hydrodynamics equations in 
one dimension. At each timestep these are first integrated by 
FCT for a half-step to evaluate time-centered fluxes; the FCT is 
then applied to a full timestep. 

Multiple dimensions are treated through directional timestep 
splitting. In multiple dimensions, the code has overall second- 
order accuracy in regions where the flow is continuous and pro- 
vides a sharp, non-oscillating solution near flow discontinuities. 
To avoid excessive temperature fluctuations at shocks, the gas 
density is smoothed over the seven nearest nodes (one cell in 
each direction) when estimating the temperature from the to- 
tal energy, velocity and density. This smoothing is done only 
for temperature estimates. Tests of the code and applications to 
cosmological problems may be found in Klypin et al. (1992) 
and Yepes ef fl/. (1995, 1996). 

The code has been fully parallelized for various shared mem- 
ory platforms. The simulation reported here was performed 
on the CRAY-YMP at CIEMAT (Spain) using 4 processors si- 
multaneously. Due to memory limitations, the supplied initial 
conditions were resampled from the original 256^ grid onto a 
coarser grid with 160 cells and particles per dimension. The 
initial particle positions were set up at z = 20 using cloud-in- 
cell interpolation of the original displacement field and veloci- 
ties were assigned by means of the Zel'dovich approximation. 
The cluster center was found iteratively from the center of mass 
of the particle distribution in spheres of radius equal to 2 cells. 

3.3. Dark matter only 

Warren - Tree 

This dark matter only simulation was carried out using a par- 
allel treecode (Warren and Salmon 1993, Warren and Salmon 
1995) on 128 processors of the 512 processor Intel Delta at 
Caltech. The algorithm computes the forces on an arbitrary 
distribution of masses in a time which scales with the particle 
number, A^, as NlogN. The accuracy of the force calculation 
is analytically bounded, and can be adjusted via a user defined 
parameter. 

Initial conditions were obtained by perturbing the masses of 
the particles in proportion to the values of the supplied initial 
density field, starting at a redshift of 63. Growing mode veloc- 
ities were assigned using the Zel'dovich approximation. The 
initial conditions were coarsened at radii exceeding 24 Mpc 
by grouping cells 8 to 1, resulting in a total of 5 340952 dark- 
matter particles in the simulation. Periodic boundary conditions 
were implemented by using the treecode to obtain forces from 
the 26 neighboring cube images, and an analytic treatment for 
the remainder. The initial portion of the simulation (to z = 9) 
was performed with a comoving Plummer softening of 50 kpc, 
and a logarithmic timestep. At z = 9, the softening was fixed at 
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5 kpc in physical coordinates, and a global timestep of .005 Gyr 
was used, resulting in 2550 total timesteps for the simulation. 
The center was defined as the particle with the highest density, 
smoothed with a spline kernel of width 20 kpc. 

The upper limit for the fractional interaction force error was 
set to 0.005. In the initial stages of the simulation, about 1200 
interactions per particle were computed. Near the end, this had 
grown to about 2300 interactions per particle. In terms of wall- 
clock time, this corresponded to about 144 seconds per timestep 
initially, and 215 seconds per timestep towards the end, repre- 
senting a sustained throughput of roughly 1.5 Gflops on the 128 
processors. 

4. RESULTS 

We first present a qualitative comparison of the results of the 
different simulations using a selection of the images. We then 
discuss quantitative results for the bulk properties and radial 
profiles of the clusters. 

4.1. Images 

The images display projections of the following quantities in 
the inner 8 Mpc cube of each simulation: 

(i) Dark matter density at z = (Figures 1 and 8) and z = 0.5 
(Figure 2), 

(ii) Gas density at z = (Figure 3) and z = 0.5 (Figure 4), 

(iii) Gas temperature at z = (Figure 5) and z = 0.5 (Figure 6), 

(iv) X-ray luminosity at z = (Figure 7). 
In Figures 1-6, the standard smoothing (250 kpc) was used, 
whereas in Figures 7 and 8 the optimal smoothing chosen by 
each simulator was used (see Table 1). The time elapsed be- 
tween the two epochs shown in the Figures is 6 x 10^ years, 
almost exactly half the dynamical time of the final cluster (de- 
fined as tdyn = 27r(r^/GM)'''^ where r^ is the virial radius and 
M the mass within it). Warren's dark matter simulation is il- 
lustrated only in Figures 1, 2 and 8, in the bottom right hand 
corner occupied in the remaining figures by Wadsley's simula- 
tion, which was the last to be completed. Wadsley's dark matter 
distribution has a very similar appearance to Couchman's. 

Dark matter density. 



All simulations show a pleasing similarity in the overall ap- 
pearance of the projected dark matter density at the final epoch 
(Figures 1 and 8). The size, shape and orientation of the main 
mass concentration are very similar in all cases. The cluster is 
elongated in the direction of a large filament - clearly visible at 
z = 0.5 (Figure 2) - along which sublcumps are accreted onto 
the cluster There are, however, noticeable differences at both 
epochs in the substructures present in the various simulations. 
These differences are due to discrepancies in the boundary con- 
ditions (assumed to be isolated in Navarro, Steinmetz and Wad- 
sley and periodic in the rest), in the treatment of tidal forces, 
and in the effective timing within the different simulations. 

The models have been evolved for at least 21 expansion fac- 
tors and inaccuracies in the initial conditions, tidal forces, or in- 
tegration errors in the linear regime lead to a lack of synchrony 
at later times. These timing discrepancies are manifest in the 
differing relative positions of some subclumps at z = 0.5 and 
are still apparent at z = 0. For example, there are two distinct 
substructures in Figures 1 and 8 to the NW of the main clump 
at z = in slightly different positions in Bryan, Cen, Jenkins, 
Owen, Pen and Warren. In Couchman, Evrard, Gnedin, Stein- 
metz and Yepes, one of these substructures is already merging 
with the central clump, while in Navarro both of them have 



merged. The differences in the overall shape and orientation 
of the main concentration are largely due to a mismatch in the 
epoch at which substructures are accreted. 

With a uniform, 250 kpc, smoothing significant noise is vis- 
ible in Owen. Figure 8 shows the dark matter distribution at 
z = 0, this time using the smoothing considered as optimal by 
each simulator to display what they considered to be real struc- 
ture. There is a larger variety of structure in these high reso- 
lution images than in the uniform smoothing case of Figure 1 . 
The simulations of Jenkins and Warren which have the largest 
number of resolutions elements also have the largest number of 
satellite structures. Varying numbers of these can be seen in 
other images, although because of the slight timing differences 
they often appear in different locations. The low and interme- 
diate resolution simulations of Cen, Gnedin, Owen and Yepes 
look quite similar when the standard smoothing or the smooth- 
ing of choice is used. 

Gas density. 

At the present epoch, the gas in all simulations (Figure 3) is 
rounder than the dark matter - a manifestation of the isotropic 
gas pressure - and has only a residual elongation along the ac- 
cretion filament. Most of the secondary clumps seen in the dark 
matter are also seen in the gas, but they are clearly more diffuse. 
At z = 0.5 (Figure 4), shortly before the final large merger, the 
timing differences discussed above are quite apparent. In some 
cases, the final major merger is already quite advanced but in 
others two large subclumps are still clearly visible. 

With this smoothing, the sampling in Owen is poor and 
Yepes' comparatively low resolution is more apparent than in 
the corresponding dark matter image. In Owen's case, the un- 
derlying asphericity of the SPH sampling is lost when a fixed 
smoothing length is used; less noisy images result when the ge- 
ometry of the hydrodynamical sampling is maintained, as in the 
X-ray image in Figure 7 below. 

Gas temperature. 

The gas temperature images (Figures 5 and 6) show the most 
interesting differences between the simulations. These are par- 
ticularly striking at z = 0.5 when the slight timing differences 
apparent in the dark matter and gas density plots produce quite 
dramatic differences in temperature structure. In particular, in 
the simulations by Jenkins, Navarro and Steinmetz, in which 
the final major merger has not yet occured at z = 0.5, an an- 
nulus of shock-heated material is evident between the two ap- 
proaching clumps, surrounding the axis of collision. In Evrard 
and Owen, the merger is further advanced, but some residue 
of the annular structure remains. The simulations of Bryan, 
Cen, Gnedin, and Yepes are yet further advanced and while 
their temperature plots show similar departures from symme- 
try, the annular structure is no longer evident. In Couchman 
and Wadsley, the merger is nearly complete and the tempera- 
ture distribution appears close to spherically symmetric. 

At z = the plots are broadly similar and most of the temper- 
ature structure, both inside the cluster and in the outer regions 
of the clusters, reproduces amongst the different simulations. 

X-ray surface brightness. 

Like the dark matter distributions in Figure 8, the X-ray 
surface brightness images in Figure 7 were generated using 
each simulator's smoothing of choice. Since the X-ray sur- 
face brightness is calculated by integrating Cx = p^T^I^, these 
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images are similar to those of the surface density (Figure 3), 
except that higher weight is given to the central parts of the 
clumps. As a result, the central intensity is strongly depen- 
dent on resolution. In the region where the bulk of the X-rays 
are produced, most images are quite similar. However, there 
are large differences in the number and brightness of satellite 
structures. Jenkins and Wadsley produced the largest number 
of such structures, followed by Bryan and Steinmetz. Again, 
the substructures appear in different places because of the tim- 
ing discrepancies discussed above. 

The high resolution grid codes of Bryan and Pen, and the in- 
termediate resolution grid code of Gnedin, show sharper struc- 
ture in the low density regions than the SPH codes, a reflec- 
tion of their better treatment of shocks in low density regions. 
Bryan's, Pen's, and all but Owen's SPH codes have high central 
resolution and generally produce brighter central regions than 
the low and intermediate resolution grid codes. 

4.2. Global properties 

All simulators were requested to calculate global properties 
of their clusters at the final time. The cluster was defined to 
be the material lying within a spherical region around the cen- 
ter, of radius such that the mean enclosed mass density is 201 
times the critical value. In practice, each simulator was left free 
to choose a preferred algorithm for locating the cluster center 
since this choice is part of the uncertainties that we are trying 
to assess. (The algorithms used in each case are described in 
Section 3.) The mean value of the cluster radius, averaged 
over all simulations, was r2oo = 2.70 Mpc, with an rms scat- 
ter of 0.04. We display the global properties of the clusters in 
graphical form in Figure 9, and discuss the different quantities 
one at a time. In each panel in this figure, the simulations are 
arranged, from left to right, in order of decreasing resolution, 
which is taken to be the maximum of the spatial resolution for 
the gas at cluster center (column 3 of Table 1) and the gravita- 
tional force resolution (column 6 of Table 1). Open circles are 
used to represent SPH simulations and filled circles grid based 
simulations. 

The simplest property of the cluster is its total mass, and all 
simulations agreed on a value just over 1O'^M0, to within better 
than 10%. Differences arise from resolution and timing effects. 
Among the grid-based codes, there is a clear trend of decreas- 
ing mass with decreasing resolution. The timing discrepancies 
discussed in Section 4. 1 affect the position of substructures and 
thus the estimates of the cluster mass. For example, among 
the SPH codes, Jenkins, Navarro, and Steinmetz find slightly 
smaller masses than the others because of a significant lump 
clearly visible in Figures 7 and 8 which falls just outside the 
cluster boundary in their simulations but just inside it in the 
others. This is a result of the late formation of their clusters, 
apparent in Figure 2. 

The velocity dispersion of the dark matter particles within the 
cluster is also reproduced to better than 10% in all simulations 
except Yepes' whose low value reflects the low force resolution 
of his simulation. (The quantity plotted is the one-dimensional 
velocity dispersion, calculated as cr/\/3, where a is the full 
three-dimensional velocity dispersion in the rest frame of the 
cluster) There is some tendency for the simulations which pro- 
duced the largest total masses also to give the largest velocity 
dispersions, but the correspondance is not exact, reflecting the 
fact that the cluster is not in virial equilibrium to better than 
10% and that the actual virial ratio depends on the detailed po- 
sitions of infalling clumps. 
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Fig. 9. Global properties of the cluster in the various simulations. All quantities are 
computed within the virial radius. From top to bottom, the left column gives the values of: 
(a) the total cluster mass; (b) the one-dimensional velocity dispersion of the dark matter; 
(c) the gas mass fraction; (d) the mass-weighted gas temperature. Also from top to bottom, 
the right column gives the values of: (a) /3 = fj.mp(Ti)^/3kT; (b) the ratio of bulk kinetic to 
thermal energy in the gas; (c) the X-ray luminosity; (d) the axial ratios of the dark matter 
(circles) and gas (crosses) distributions. In each panel, the models are arranged, from left 
to right, in order of decreasing resolution which is taken to be the maximum of the spatial 
resolution for the gas at cluster center (column 3 of Table 1) and the gravitational softening 
length (column 6 of Table 1 ). Open circles are used to represent SPH simulations and filled 
circles grid based simulations. 

As might be expected, although still quite good, the agree- 
ment on the properties of the gas is noticeably worse. The total 
amount of gas is most interestingly expressed as a fraction of 
the total cluster mass. The overall gas fraction in the model uni- 
verse is 10% and the highest resolution grid simulations find a 
cluster gas fraction which is almost exactly equal to this. Lower 
resolution grid simulations find progressively smaller gas frac- 
tions, reflecting differential resolution effects in the treatment 
of dark matter and gas in these codes. There is excellent overall 
agreement among the SPH models except for the lowest resolu- 
tion one: everyone except Owen finds a gas fraction very close 
to 0.09. There seems to be a systematic offset between the SPH 
models and the highest resolution grid models but there is no 
clear indication of what may be causing this difference. We ex- 
plore this issue further in the next section by considering the 
radial dependence of the gas fraction. 

As a measure of temperature, all simulators calculated a 
mean, mass-weighted temperature for all the gas within the 
cluster Everyone found a value between 4.9 and 5.9 x lO^K. 
Figures 5 and 6 suggest that some of the differences result from 
the timing differences discussed above which produce slightly 
different histories for the clusters just prior to z = 0. It is encour- 
aging that the rms scatter in the measured mean temperature is 
less than ±7%. An even smaller scatter, ±5%, is found for the 
ratio of specific dark matter kinetic energy to gas thermal en- 
ergy, (3 = ^nipaj^f^/'ikT, if Yepes' low value is excluded. (Here 
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/i denotes the mean molecular weight and nip the proton mass.) 
This is further evidence that the differences in velocity disper- 
sion and temperature result from slightly different dynamical 
histories rather than from differences in the treatment of the gas. 
Most simulations give /? ~ 1.17, indicating that non-thermal or 
bulk turbulent motions contribute to the support of the gas. The 
ratio of bulk kinetic to thermal energy in the gas, plotted in the 
next panel of Figure 9, is indeed about 15%. This ratio corre- 
lates well with j3, but with small residuals that reflect slight de- 
partures from virial equilibrium. The agreement on the values 
of P and Ukin/Utherm indicates that the shock capturing proper- 
ties and the efficiency with which infall kinetic energy is ther- 
malized in shocks is similar in the SPH and grid-based codes, at 
least in the regime explored in this simulation. Most simulators 
obtained values of about 3.6 x lO^^ergs for the bulk kinetic en- 
ergy, except Gnedin and Wadsley who found values about 25% 
smaller. Lower turbulent energies probably result from some 
combination of smaller noise-induced motion, greater viscous 
damping, and, in Wadsley's case a more dynamically advanced 
state, but the actual factors responsible in each case are unclear. 

There is substantially less consensus about the estimated 
X-ray luminosities of the cluster, calculated approximately as 
J p^T'^'^dV, and so given in units of MqK'/^Mpc"-'. The values 
found span a range of a factor of ~ 10, or a factor of ^ 5, if we 
exclude Yepes' low resolution model. Resolution effects also 
account for the small values obtained by Cen and Owen. The 
largest values were obtained in the intermediate and high res- 
olution grid-based simulations of Gnedin and Pen respectively, 
and in Evrard's SPH simulation. Evrard's and Wadsley's mod- 
els produced larger X-ray luminosities than other SPH simula- 
tions because their clusters are in a slightly more advanced (and 
more active) dynamical state. The higher luminosities from the 
high resolution grid-based codes are due to their slightly more 
concentrated central gas distributions (see Figure 3 and Sec- 
tion 4.3 below). Because the total X-ray luminosity is sensitive 
to the structure of the inner few hundred kiloparsecs of the clus- 
ter, it fluctuates quite strongly in time and is very sensitive to 
simulation technique. 

As a final test, we compare the shapes of the clusters, as 
measured by the inertia tensors, X = ^ m,x,x,/ ^ m,-, for both 
the dark matter and the gas within the spherical region which 
defines the cluster. We label the eigenvalues of this tensor 
a^ > b^ > c^, and define the axial ratios to be a/b and a/c. As 
can be seen in Figures 1 and 3, the cluster is aspherical, with the 
orientation of its longest axis reflecting its formation by infall 
along a filament. The axial ratios shown in Figure 9 show, in 
fact, that the cluster is triaxial. The dark matter distribution (cir- 
cles in Figure 9) is considerably more aspherical than the gas 
(crosses) in all cases except Gnedin's. There is generally good 
agreement amongst the different simulations, although there are 
a few anomalies. For example, Owen finds the smallest axial 
ratios for the dark matter distribution even though the inner re- 
gions of his cluster appear quite elongated in Figure 1 . This is 
probably because of the relatively large contribution to T from 
the largest infalling clump which lies close to raoo in his simu- 
lation. Gnedin's dark matter and gas distributions are consider- 
ably more aspherical than the rest. 

4.3. Radial profiles 

In order to perform more detailed quantitative comparisons 
of cluster structure, each simulator was asked to provide the 
radial profiles of a number of cluster properties. These were 
averaged in a specified set of spherical shells centered at the 



position deemed by each simulator to be the cluster center (cf 
§2.2). Simulators were also asked to specify the effective res- 
olution of their simulation, and throughout this section we plot 
data for each model only at radii larger than this. Most simu- 
lators identified the resolution of their simulation with the ef- 
fective gravitational force resolution (see Table 1), but the fol- 
lowing specified different values: Cen (200 kpc), Owen (500 
kpc). Pen (50 kpc), Steinmetz (35 kpc), and Yepes (400 kpc). 
In the figures that follow, the data points are slightly displaced 
from the bin centers for clarity, and we use a solid line to show 
the mean profile obtained by averaging all the data plotted in 
each bin. At the top of each diagram we plot the residuals from 
this mean, defined in most cases as Iwc- < Inx >, where x is the 
property of interest and the brackets denote the average of the 
data points in each bin. This definition applies to the residuals 
of all the radial profiles, except those of the normalized baryon 
fraction (Figure 13), which we define as (x- < x >)/ < x >, 
and those of the radial velocity profiles (Figures 14 and 15) and 
the entropy (Figure 18) which we define as x- < x >. The sam- 
pling errors in the estimates of the different cluster properties 
are largest in the innermost bin plotted where they are typically 
less than about 10%. 
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Fig. 10. Projected dark matter density at z = 0. The images, covering tlie inner 8 Mpc of 
eacli simulation cube, have been smoothed using the standard Gaussian filter of 250 kpc 
half-width described in the text. Wadsley's simulation, not shown here or in Figure 2, has 
a similar appearance to Couchman's. 

We begin by comparing the dark matter density profiles plot- 
ted in Figure 10. In general there is very good agreement be- 
tween the different calculations over the regions resolved by 
each simulation. The two highest resolution models, Warren's 
pure N-body simulation, and Jenkins' AP^M/SPH simulation, 
agree extremely well at all radii. The residuals plot shows that 
all simulations agree to within ±20% at all radii. The dark mat- 
ter profile in this cluster is well fit by the analytic form proposed 
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by Navarro, Frenk and White (1995), all the way from 10 kpc 
to 10 Mpc. This fit is shown as a dashed line in Figure 10 and 
corresponds to a value of the concentration parameter, c = 7.5, 
appropriate to a typical isolated halo of this mass in an il = 1 
CDM model (Navarro, Frenk & White 1997). Only in the very 
center is there a slight indication that the true profile, as defined 
by Warren's model, might be steeper than the analytic form. 

The dark matter velocity dispersion profiles in Figure 1 1 con- 
firm that all the codes give very similar results for the dynamical 
properties of the dark matter. (Note the very different dynamic 
ranges in Figures 10 and 11.) Except for the last bin, which 
is particularly affected by noise arising from subclustering, the 
scatter in the velocity profiles is comparable to that in the mass 
profiles, about 20%. The velocity dispersion profile rises near 
the center, has a broad peak around 100-500 kpc, and declines 
in the outer parts. Warren's N-body model and all the high res- 
olution SPH and grid models resolve the inner rising part of the 
velocity dispersion profile. 
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Fig. 11. The dark matter velocity dispersion profile. The quantity shown in the main 
plot is the one-dimensional velocity dispersion, calculated as (t/v3, where a is the full 
three-dimensional velocity dispersion. See the caption to Figure 10 for a description of the 
symbols and other details of the plot. 

The agreement of the gas density profiles (Figure 12) is 
less good but is still quite impressive. In the fixed Eulerian 
grid models of Cen and Yepes, the gas density at their inner- 
most point is somewhat low, whereas the corresponding dark 
matter matter densities agree well with the other calculations. 
This shows, unsurprisingly, that the resolution of such codes is 
somewhat poorer for the hydrodynamics than for the N-body 
dynamics. Bryan's multilevel grid code produces results that 
agree quite well with other high resolution models, except that 
his two innermost points lie slightly below those of the highest 
resolution SPH models. Gnedin's and Pen's variable resolu- 



tion grid codes give higher than average gas densities in the 
200 to 600 kpc range and, as a result, their density profiles 
are steeper than the others. Pen's simulation produced a more 
pronounced core structure within about 150 kpc than all other 
models while Gnedin's simulation shows the largest departures 
from the mean profile. There is also significant scatter among 
the results of the various SPH codes, Evrard finding system- 
atically high gas densities in the inner cluster and Couchman 
finding systematically low values. It seems likely that at least 
some of these differences result from the differences in the tim- 
ing of cluster collapse noted in Section 4.1, rather than from 
differences in the treatment of the hydrodynamics between the 
different techniques and implementations although the two can- 
not be clearly disentagled. For example, Couchman's gas den- 
sities between 80 and 200 kpc are somewhat lower than those 
obtained by Jenkins using the same SPH implementation but 
different numbers of particles and integration parameters. The 
dashed line in Figure 12 is the mean dark matter density profile 
reproduced from Figure 10. It is interesting that this profile is 
substantially steeper than the mean gas density profile in the in- 
ner cluster, indicating that the gas has developed a much more 
clearly defined "core" than the dark matter. 
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Fig. 12. The gas density profile. The dashed line shows the mean dark matter density 
profile reproduced from Figure 10. See the caption to Figure 10 for further details. 

The discrepancies between the various dark matter and gas 
density profiles become clearer when we examine the variation 
of the normalized gas fraction, T = M^as{< '')/i^bMtot(< r)), 
with radius (Figure 13). There is considerable scatter (~ 50%) 
in this gas fraction over the 100 kpc to 1 Mpc range, Gnedin 
finding the highest values and Couchman, Cen, and Yepes the 
lowest. Well inside the virial radius, three of the grid models, 
Bryan's, Pen's, and especially Gnedin's, rise above T = 1, the 
mean for the simulation as a whole. At larger radii. Pen's val- 
ues fall slightly below this mean, while Bryan's and Gnedin's 
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remain slightly above unity well beyond the virial radius. By 
contrast, none of the SPH models ever rise above T = 1 and 
they all give very similar results at the virial radius, similar also 
to the Eulerian models of Cen and Yepes. Although relatively 
small, these discrepancies appear to reflect a systematic differ- 
ence in the final distribution of the gas between the three inter- 
mediate and high resolution grid simulations (Gnedin, Pen and 
Bryan) and the rest of the models. Particularly puzzling is the 
excursion towards large values of T seen by Gnedin at r ~ 1 
Mpc and the fact that Bryan and Gnedin obtain values of T > 1 
beyond three virial radii where one might expect the mixture of 
dark matter and gas to attain the mean universal value. Note, 
however, that the deviations from T = 1 at large radii are quite 
small. 

The infall patterns of dark matter and gas around the clus- 
ter (i.e. the run of peculiar radial velocities) are illustrated in 
Figures 14 and 15. For the dark matter, the radial velocity pro- 
files are quite similar: net infall is seen in most models beyond 
~ 700 kpc, except in Pen's case, in which the radial velocity 
remains close to zero until about twice this radius. There are 
larger differences in the radial velocity profiles for the gas. In 
most models the gas is infalling over the same range of radii as 
the dark matter, but at a slightly lower speed. Gnedin's model 
is anomalous in this respect. In Pen's simulation, on the other 
hand, there is a small net outflow of gas at ^ 1 Mpc. The scat- 
ter in the radial velocity profiles beyond ^ 1 Mpc is about 200 
km s~' . 
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Fig. 13. The radial dependence of the gas fraction. The quantity plotted is the gas 
fraction normalized to the value for the simulation as a whole (10%). See the caption 
to Figure 10 for further details, but note that in this Figure the residuals are defined as 
(T- < T »/ < T >. 



Fig. 14. The radial velocity profile of the dark matter. Velocities are computed in the rest 
frame of the cluster and do not include the Hubble expansion. See the caption to Figure 10 
for further details, but note that in this Figure the residuals are defined as v^— < v^ >. 
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Fig. IS.The radial velocity profile of the gas. Velocities are computed in the rest frame 
of the cluster and do not include the Hubble expansion. The dashed line is the dark matter 
radial velocity profile from Figure 14. See the caption to Figure 10 for further details, but 
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note that in this Figure the residuals are defined as iv— < vv >■ 

The differences between the codes become most obvious if 
we look at various thermodynamic properties of the gas. Fig- 
ures 16 through 19 show radial profiles for the pressure, tem- 
perature, entropy and X-ray emissivity. (Of these, only the 
temperature was calculated directly by the simulators; the other 
quantities were derived from the binned values of temperature, 
density, and X-ray luminosity.) Agreement among the pres- 
sure profiles is reasonably good. This reflects the good agree- 
ment of the dark matter density profiles and so of the gravita- 
tional potential wells, together with the fact that the intracluster 
gas is close to hydrostatic equilibrium in all cases. In the cen- 
tral regions, simulations which give flatter gas density profiles 
(e.g. Bryan, Pen, Couchman, and Cen) can be clearly seen to 
give correspondingly higher temperatures, as required in order 
to maintain comparable pressure gradients. The differences in 
the temperature profiles appear substantially larger than in other 
quantities, but this is in part a reflection of the smaller dynamic 
range of the plot. There is a suggestion that the temperature 
structure in the inner parts may be systematically different in 
the SPH and grid models. In the former, the temperature profile 
is flat or slowly declining towards the inner regions, but in the 
grid simulations, the temperature is still rising at the innermost 
point plotted, a trend that is particularly noticeable in Bryan's 
simulation. In the outer cluster, the temperature profile drops 
substantially in all cases. 
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Fig. 16. The gas pressure profile. See the caption to Figure 10 for further details. 
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Fig. 17.The gas temperature profile. The quantity plotted is the mass-weighted tempera- 
ture. See the caption to Figure 10 for further details. 

The entropy and X-ray luminosities show the patterns ex- 
pected for quantities derived in a simple way from the density 
and temperature of the gas. Note, in particular, that the entropy 

2/3 

(defined as s = ln(T /pgas) decreases systematically towards the 
center in all the SPH models, but that this decline is less obvious 
in the grid models and is, in fact, absent in the central parts of 
Bryan's simulation in which the entropy remains approximately 
constant within ^ 200 kpc. This difference might reflect differ- 
ences in the way in which shocks are treated in the SPH and grid 
codes; however, the effect is small and occurs at the resolution 
limit of the grid simulations. Finally, Figure 19 shows the con- 
tribution to the total X-ray luminosity per logarithmic interval 
in radius, Airr^Cx- Most of the X-ray luminosity is produced 
in the radial range 200-500 kpc. This region was well resolved 
by all the SPH simulations and by Bryan's model which agrees 
remarkably well with them. On the other hand, Cen, Owen and 
Yepes did not resolve this radial range and, as a result, their to- 
tal luminosities ended up being smaller than average. The large 
gas densities found by Gnedin and Pen in this range account for 
their larger than average total luminosities (see Figure 9). Thus, 
the large scatter in total X-ray luminosity seen in Figure 9 re- 
sults partly from the low resolution of some of the models and 
partly from the unusually high gas densities found by the two 
deformable grid models. Between ^ 1 Mpc and the virial radius 
the different models agree better although the scatter in AirPCx 
is still larger than in all other properties. 
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Fig. IS.The radial variation of the gas entropy. The entropy is defined as s = ln(T /p'^ 
See the caption to Figure 10 for further details. 
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Fig. 19. The X-ray luminosity profile. The quantity plotted is 47vr^Cx, where Cx de- 
notes the X-ray luminosity density in each bin. See the caption to Figure 10 for fuither 
details. 



5. DISCUSSION AND CONCLUSIONS 

We have simulated the formation of an X-ray cluster in a 
cold dark matter universe using 12 different cosmological gas 
dynamics codes that span the range of numerical techniques 
and implementations presently in use. This comparison aims to 
assess the reliability of current cosmological simulations in the 
regime relevant to the bulk of the gas in galaxy clusters, and to 
set a standard against which future techniques may be tested. 
(T he initial conditions and a selection of r esults are available 
at http://star-www.durac.uk^csf/clusdata/ or by request from 



CSF.) Because our goal is to compare results in a realistic sit- 
uation, only the initial conditions were specified. Other vari- 
ables such as resolution, the treatment of boundary conditions, 
and integration parameters were left to the discretion of simula- 
tors. Our comparison therefore encompasses not only the cur- 
rent variety of hydrodynamics techniques, but also the different 
choices commonly made by individual authors for these vari- 
ables. Seven of the codes used for this comparison are based 
on the SPH technique, while the other five are based on grid 
methods, employing either a fixed, a deformable or a multi- 
level mesh. The resolution of the simulations varied from a few 
tens to a few hundreds of kiloparsecs. Although there is, of 
course, no guarantee that any of the calculations gives the cor- 
rect solution to the problem, the agreement among the various 
simulations was better than a pessimist might have predicted. 
Nevertheless, some important differences do exist. 

The simulated cluster was chosen to have a mass comparable 
to the Coma cluster at the present day and to appear fairly re- 
laxed on visual inspection. We assumed cold dark matter initial 
conditions, fl=l, h = 0.5, and a global baryon fraction of 10%. 
The cluster formed by the merging of subclumps infalling along 
a filament and experienced a final major merger between z = 0.5 
and z = 0. The final cluster mass, virial radius, gas fraction 
within the virial radius, and mass-weighted gas temperature, 
averaged over all the simulations, and the standard deviation 
in each quantity areM= 1.1 x lO'^M©, ctm = 0.05 x lO'^M©; 
r, = 2.70 Mpc, an. = 0.04 Mpc; fi, = 0.092, cr/„ = 0.006; and 
r = 5.4 X lO^K, (Tt = 0.34 x lO^K respectively. 

The properties of the cluster dark matter are gratifyingly sim- 
ilar in all the models. The total mass and velocity dispersion 
agree to better than 5%. The dark matter density and veloc- 
ity dispersion profiles are also similar and match the result of 
a higher resolution dark matter only simulation. Over the re- 
gions adequately resolved in each simulation, the scatter in 
these quantities, relative to the mean profile, is less than about 
±20% per logarithmic bin in radius. For the most part, this 
scatter seems to be due to a slight asynchrony in the evolution- 
ary state of the models introduced by inaccuracies in the initial 
conditions, the treatment of boundary conditions and of tidal 
forces and the integration procedure. Thus, small subclumps 
often appear at different positions and the timing of the final 
major merger differs slightly in the different models. 

There is less agreement on the gas properties of the cluster, 
although in most cases they are quite similar For example, all 
models agree to 10% (rms) on the gas mass and baryon frac- 
tion within the virial radius. In all the SPH and all but one of 
the grid models, the gas is slightly more extended than the dark 
matter. The scatter relative to the mean in logarithmic radial 
bins seldom exceeds 20% in the case of the density profile and 
30% in the case of the gas mass fraction. There is no obvious 
systematic difference in the gas density profiles produced by 
the SPH and fixed grid models, but the deformable grid models 
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produced somewhat larger core radii. At large distances from 
the cluster, some of the grid based models rise slightly above 
the theoretical expectation of a universal baryon fraction. 

The mean (mass-weighted) gas temperature is reproduced to 
within 6% (rms) by all the codes. The ratio of specific dark mat- 
ter kinetic energy to gas thermal energy is reproduced to a simi- 
lar accuracy. The scatter per logarithmic bin in the temperature 
profiles falls in the ±20% band and is only slightly larger for the 
pressure profile. In the central regions (r < lOOkpc), however, 
the SPH codes produce a flat or slightly declining temperature 
profile while all the grid codes produce a temperature profile 
that is still rising at the resolution limit. The entropy of the gas 
declines continuously from the virial radius to the resolution 
limit but there is a suggestion that the entropy in the grid codes 
may bottom out at small radii while it continues to decrease in 
the SPH codes. 

Amongst all the properties we have examined, the largest dis- 
crepancies occur in the predicted cluster X-ray luminosity. This 
is proportional to the square of the gas density and so is strongly 
dependent on resolution and is further affected by variations in 
the potential produced by the exact timing of the final major 
merger. The large range of effective resolution in the various 
models gives rise to a factor of 10 variation in the total X-ray 
luminosity. The luminosity per logarithmic interval in radius 
peaks just outside the gas core radius. The eight simulations 
that resolved this region (all the SPH and two of the grid mod- 
els) show a much narrower spread in total X-ray luminosity, 
amounting to a factor of 2.6 (or 1.8 if the most extreme model 
is excluded.) 

We conclude that the different approaches to modelling 
shocks and other hydrodynamical processes implicit in the di- 
verse techniques employed in this comparison give, in most 
cases, fairly consistent results for the dynamical and thermo- 
dynamical properties of X-ray clusters. Variations introduced 
by differences in the internal timing of the simulations tend to 
be at least as important as variations in the treatment of hy- 
drodynamics. An illustration of this is the comparison of the 
simulations of Couchman and Jenkins who used serial and par- 
allel versions of essentially the same code, but with different 
numbers of particles and other simulation parameters. The fi- 
nal temperatures of their clusters differed by 15%, the X-ray 
luminosities by 20%, and the bulk gas kinetic energy by 15%. 

The conclusions discussed in this paper apply exclusively to 
the particular case of a non-radiative gas. They cannot be ex- 
trapolated to other regimes such as that appropriate to galaxy 
formation where gas cooling is a dominant process. The behav- 
ior of the gas in this situation is determined by the resolution 
of the simulation because the cooling rate depends strongly on 
gas density. In addition, simulations of galaxy formation often 
include algorithms to convert cold gas into stars and to model 
the feedback processes associated with star formation. A com- 
parison of galaxy formation simulations, analogous to the com- 
parison of X-ray cluster formation carried out in this paper is 
clearly desirable, but would be much more complex to imple- 
ment in practice. The level of agreement among the techniques 
currently in use for studies of galaxy formation remains an open 
question. Perhaps a simpler, but instructive, next step might be 
a comparison of simulations of high redshift gas clouds (the 
"Lyman-a forest" clouds), in which gas cooling is less impor- 
tant. 

Based on the overall consistency of the simulations discussed 
in this paper, we can draw a number of conclusions regarding 
the properties of the simulated cluster which we expect to be 



typical of near equilibrium, massive clusters in a CDM uni- 
verse. Some of these conclusions mirror those found in ear- 
lier work (eg. Evrard 1990, White et al. 1993, Bryan et al. 
1994, Cen & Ostriker 1994, Kang et al. 1994, Navarro, Frenk 
& White 1995, Anninos & Norman 1996, Bartelman & Stein- 
metz 1996.) 

1) The dark matter distribution in our simulated cluster is elon- 
gated along the direction of the dominant large filament along 
which subclumps were accreted onto the cluster. The final gas 
distribution is rounder than the dark matter distribution and, as 
a result, the direction of the filament is difficult to identify in an 
X-ray image. 

2) The dark matter density profile in the cluster is well fit by the 
analytic form proposed by Navarro, Frenk and White (1995), 
all the way from 10 kpc to 10 Mpc. This form has a radial de- 
pendence close to r~' in the inner regions, steepens to r~^ at 
intermediate radii and falls off like r~^ in the outer parts. The 
corresponding velocity dispersion profile rises from the center 
outwards, has a broad maxium around 100-500 kpc, and de- 
clines in the outer parts. 

3) Although the gas is close to hydrostatic equilibrium through- 
out the cluster, mergers disturb both the gravitational potential 
and the dynamical state of the gas. Bulk motions make a small 
but significant contribution to the support of the gas: the kinetic 
energy in bulk gas motions is about 15% of the thermal energy 
of the gas. The quantity /? = jinipCrj^f^/ikT has a mean (aver- 
aged over all the simulations except the lowest resolution one) 
of 1.15 with an rms scatter of only 0.05. 

4) The radial density profile of the gas in the highest resolu- 
tion simulations develops a "core radius", ie. a region in which 
the slope of the profile flattens rapidly. This change of slope oc- 
curs at r<250 kpc and, inside this radius, the gas density profile 
is significantly flatter than the dark matter density profile. The 
gas fraction within a given radius therefore rises from the center 
outwards and, at the virial radius, the mean over all the simula- 
tions is 0.92 of the global value, with a fractional rms scatter of 
only ±0.065. The temperature distribution in the inner parts has 
an approximately flat profile and, beyond a few hundred kpc, it 
begins to decline so that at the virial radius, the temperature is 
^ 30% of the central value. 

5) A reliable estimate of the cluster X-ray luminosity requires 
resolving the radial range 200-500 kpc, or 5%-20% of the virial 
radius, where the X-ray luminosity per logarithmic interval in 
radius peaks. Even when this is possible, the strong sensitivity 
of X-ray luminosity to local variations in gas density leads to 
a spread of a factor of ~ 2 in the predicted X-ray luminosity. 
These variations are due to a variety of numerical effects, and 
a factor of 2 uncertainty is a realistic estimate of the accuracy 
with which cluster X-ray luminosities can be predicted with the 
present generation of techniques and computing resources. 
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Fig. 1 . — Projected dark matter density at z = 0. The images, covering the inner 8 Mpc of each simulation cube, have been smoothed using the standard Gaussian 
filter of 250 kpc half-width described in the text. Wadsley's simulation, not shown here or in Figure 2, has a similar appearance to Couchman's. 



Fig. 2. — Projected dark matter density at z = 0.5. The images, covering the inner 8 Mpc of each simulation cube, have been smoothed 
using the standard Gaussian filter of 250 kpc half- width described in the text. 



Fig. 3. — Projected gas density at z = 0. The images, covering the inner 8 Mpc of each simulation cube, have been smoothed using 
the standard Gaussian filter of 250 kpc half- width described in the text. 



Fig. 4. — Projected gas density at z = 0.5. The images, covering the inner 8 Mpc of each simulation cube, have been smoothed using 
the standard Gaussian filter of 250 kpc half- width described in the text. 



Fig. 5. — Integrated, mass-weighted gas temperature at z = 0. The images, covering the inner 8 Mpc of each simulation cube, have 
been smoothed using the standard Gaussian filter of 250 kpc half- width described in the text. (The roughly circular "cut-out" regions 
seen in the outer parts of this and the next figure are associated with cool, infalling clumps of size comparable to the resolution of the 
smoothed image (see Figures 3 and 4); the edges are enhanced by the choice of color table.) 



Fig. 6. — Integrated, mass-weighted gas temperature at z = 0.5. The images, covering the inner 8 Mpc of each simulation cube, have 
been smoothed using the standard Gaussian filter of 250 kpc half-width described in the text. 



Fig. 7. — Projected X-ray luminosity at z = 0. The images, covering the inner 8 Mpc of each simulation cube, have been smoothed 
with the filter chosen by each author to best portray the results of each simulation (see Table 1). 



Fig. 8. — Projected dark matter density at z = 0. The images, covering the inner 8 Mpc of each simulation cube, have been smoothed 
with the filter chosen by each author to best portray the results of each simulation (ee Table 1). 
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